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Abstract. We formulate an atomistic-to-continuum coupling method based on blending atomistic 
and continuum forces. Our precise choice of blending mechanism is informed by theoretical predic- 
tions. We present a range of numerical experiments studying the accuracy of the scheme, focusing 
in particular on its stability. These experiments confirm and extend the theoretical predictions, 
and demonstrate a superior accuracy of B-QCF over energy-based blending schemes. 
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1. Introduction 

Atomistic-to-continuum coupling methods (a/c methods) have been proposed to increase the 
computational efficiency of atomistic computations involving the interaction between local crystal 
defects with long-range elastic fields [6|[7|[l5|[T9|[23|[30||3l|[42]; see [27j for a recent review of a/c 
coupling methods and their numerical analysis. Energy-based methods in this class, such as the 
quasicontinuum model (denoted QCE [43| ) , exhibit spurious interfacial forces ("ghost forces") even 
under uniform strain |8y41j . The effect of the ghost force on the error in computing the deformation 
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and the lattice stability by the QCE approximation has been analyzed in 
ment of more accurate energy-based a/c methods is an ongoing process [5 

An alternative approach to a/c coupling is the force-based quasicontinuum (QCF) approxima- 
tion [7, 11 12ti26j^0j, but the non-conservative and indefinite equilibrium equations make the iter- 
ative solution and the determination of lattice stability more challenging 12-13. Indeed, it is an 
open problem whether the (sharp-interface) QCF method is stable in dimension greater than one. 
Although some recent results in this direction exist [251, it is still unclear to what extent they can 
be extended for general atomistic domains and in the presence of defects. 

Many blended a/c coupling methods have been proposed in the literature, e.g., [l}|4| [T6||24[[37| 
38,44^. In [22^, we formulated a blended force-based quasicontinuum (B-QCF) method, similar to 
the method proposed in [26] , which smoothly blends the forces of the atomistic and continuum 
model instead of the sharp transition in the QCF method. Under the simplifying assumption 
that deformation is homogeneous, we established sharp conditions under which a linearized B- 
QCF operator is positive definite, which effectively guarantees stability of the numerical scheme. 
The required blending width to ensure positive definiteness of the linearized B-QCF operator is 
surprisingly small. In the present paper, we present focused numerical experiments to confirm 
and extend the theoretical conclusions in 20 , 22 , and in addition provide accuracy benchmarks 



similar to those in [28] . Our numerical benchmarks demonstrate that the B-QCF scheme is a 
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practical a/c coupling mechanism with performance (accuracy versus computational cost) superior 
to energy-based blending schemes. 

1.1. Summary. In section [2| we introduce the B-QCF model for a ID atomistic chain. We state 



the asymptotically optimal condition on the blending size in Theorem 2.1 and apply a uniform 



expansion to the atomistic chain in subsection 2.2, The critical strain errors between the atomistic 



and B-QCF models with different blending size are computed in this subsection. The numerical 
results perfectly match the analytic prediction, that is, the errors decay polynomially in terms of 
the blending size. 

In section [Sj we establish the B-QCF model for a 2D hexagonal lattice. We state sufficient and 
necessary conditions on the blending width under which the B-QCF operator is positive definite. 
To investigate the positive-definiteness of the B-QCF operators in 2D, we apply three different 
classes of deformations to the perfect lattice, which are the uniform expansion, two types of shear 
deformation, and a general class of homogeneous deformations. The results of 2D uniform expansion 
are similar to those of the ID example, and they agree with the theoretical conclusions well. 

The stability regions of the different models under homogeneous deformations are consistent 
with the analytic prediction. By using a small blending region, the 2D B-QCF operator becomes 
almost as stable as the atomistic model, compared to the fact that the stability region of the force- 
based quasicontinuum (QCF) method, i.e., the B-QCF method without blending region, is a proper 



subset of the fully atomistic model 12-14 . However, the stability error under shear deformation 
for the B-QCF operator seems to only depend linearly on the system size, which is observed from 
the numerical experiments. 

In section |4j we implement the B-QCF method from a practical point of view. We briefly review 
the accuracy results in terms of computational cost, i.e., the total number of degrees of freedom DoF, 
and then include some numerical experiments for a di- vacancy and microcrack to demonstrate the 
superior accuracy of B-QCF over other a/c coupling schemes that we have investigated previously 
in (281. 



2. The B-QCF Operator in ID. 

2.1. Notation. We denote the scaled reference lattice by eZ := {ei : i G Z}. We apply a macro- 
scopic strain F > to the lattice, which yields 

YF :=FeZ = {F€i)eez. 
The space U of 2A^-periodic zero mean displacements u = (u^)^gz from yp is given by 

U := <u: ue+2N = ue for £ £ Z, and J2iL-N+i Ui = 0>, 

and we thus admit deformations y from the space 

yp := {y : y = Yf + u for some u e U}. 

We set e = 1/A^ throughout so that the reference length of the computational cell remains fixed. 
We define the discrete differentiation operator, Du, on periodic displacements by 

[L)u)i := , — oo < t < oo. 

e 
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We note that {Du)^ is also 2A'^-periodic in i and satisfies the zero mean condition. We will often 
denote (-Du)^ by Du£. We then define (-D^^'u) „ and (D^^'u) „ for — oo < ^ < cxo by 



e \ / i e 

To make the formulas more concise, we sometimes denote Dui by u'^, D^'^^ui by u", etc., when there 
is no confusion in the expressions. 

For a displacement \i £ U and its discrete derivatives, we employ the weighted discrete £? and 
i"^ norms by 

/ N \ Vp 

IIuILp := e > \ufF for 1 < p < oo, llulUoo := max luA, 

\ e=-N+i / — 

and the weighted inner product for i'^ is 

N 
(u, w):= y^ £U£W£. 
e=-N+i 

2.2. The B-QCF Operator. We consider a one-dimensional (ID) atomistic chain with periodicity 
2N, denoted y G yp, under second- neighbor pair interaction. The total atomistic energy per period 
of y is given by S^iy) - e YleL-N+i fiV^i where 

N 

r{y) = e Y. [<Piy'i) + Hyi + y'i-i)] (2.1) 

e=-N+i 
for external forces /^ and a scaled Lennard-Jones type potential 18,33 cj) ^ C^(0,+oo). Implicitly 



we also assume that (t>{r), (p'{r) and (p"{r) decay rapidly as r increases, so that we only have to take 
into account first and second neighbors. 

The equilibrium equations are given by the force balance at each atom: F^ + fi = where 

The linearized equilibrium equations about yp are 

(L^u^), = /,, for i = -N + l,...,N, 
where (L^v) for a displacement v £U is given by 

[L v)^ := <f>p "2 + 4>2F ^2 • 

Here and throughout we use the notation (pp := (j)" {F) and (j)2p '■= 4>"{2F), where (f) is the potential 



in (2.1 ). We assume that (j)p > 0, which holds for typical pair potentials such as the Lennard-Jones 
potential under physically relevant deformations. Appropriate extensions of the stability results 
in this paper can likely be obtained for more general smooth deformations by utilizing the more 



technical formalism developed, for example, in 19,34,35]. 



The local QC (or Cauchy-Born) approximation (QCL) uses the Cauchy-Born extrapolation rule 



42,43 , that is, approximating y'^ + y'^_^ in (2.1) by 2y'^ in our context. Thus, the QCL energy is 
given by 

N 

£'''\y) = ^ E [Hyi) + HH)]- (2.3) 

e=-N+i 
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Then the local continuum forces F^'^^(y) are 

We can similarly obtain the linearized QCL equilibrium equations about the uniform deformation 

^^qcl^qcl^ =/f for i=-N+l,...,N, 

where the expression of (L^'^ v) „ with v G Z// is 

The blended QCF (B-QCF) operator is obtained through smooth blending of the atomistic and 
local QC models. Let /3 : M — t- M be a "smooth" and 2-periodic blending function, then we define 

i^^''(y) := PlFliy) + (1 - Pi)Ffiy), 

where f3i := (3{ei). Linearization about yp yields the linearized B-QCF operator 

(L^-^^V), := Pe{L^v)e + (1 - /3,)(L^=V),. 

Next, we define the blending width K: 

I:= {i(z {_Ar + 1, . . . AT} : < /S^^j < 1 for some j € {±1, ±2}} and K := ^I, (2.4) 

so that D^^^Pi = for ah i G {-N + 1, . . . N} \I and j G {1, 2, 3}. Thus K is the size of the 
compact support of D'^^' (3i. It is obvious that K < 2N . 

2.3. Positive-Definiteness of the B-QCF Operator. We proved in [22] that the blending 
function j5 can be chosen as a quintic polynomial such that 
(i) The jth derivatives of (3 satisfy 

p(^')/3||,^ < CpiKey^, for j = 1, 2, 3. (2.5) 

(ii) This estimate is sharp in sense that, if /3£ attains both the values and 1, then 

P^^Vlk- > {Ke)~^, for j = 1, 2, 3. (2.6) 

A linearized operator U" with w G {a, c,bqcf}, is said to be positive definite in the H^ norm or 
coercive if there exists a constant 7 > such that 

(L™u,u) >7||Du||22 VuG^/. (2.7) 

We have proved an asymptotically optimal stability condition on the blending region size of the ID 



B-QCF operator in 22 



Theorem 2.1. Let X and K he defined as in (2.4), and suppose that j3 is chosen to satisfy the 



upper bound ( |2.5[ ). Then there exists a constant Ci = Ci(C^), such that 

(L'^^^fu,u) > (co - Ci|0'2V|[i^~'/'iV'/'])||I?u||| Vu G^/, (2.8) 

where cq = min(0^, (pp + ^4'2p) is the atomistic stability constant. 

Moreover, if Pi takes both the values and 1, then there exist constants 02,0^ > 0, independent 
of I, N, (f)'p and (/>2^, such that 

(L^q^^u, u) < (co + {^2 - C3 [i^-5/2^1/2] } \cp'ip\^ \\Du\\\ Vu G U. (2.9) 
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From the conclusion of Theorem 2.1 , we can immediately get the following necessary and sufficient 
conditions on the blending width K for the operator L^^'^^ to be coercive. 

Corollary 2.1. Suppose that L^ is positive-definite and that the blending function is sufficiently 
smooth. If the blending size K satisfies K 3> N^'^ , then the B-QCF operator L^'^'^^ is positive- 
definite and this estimate is asymptotically optimal. 

2.4. ID Uniform Expansion Experiments. We conduct numerical experiments in order to 
verify our theoretical findings. More precisely, we compare the decay rates of the error in critical 
strain as computed by B-QCF with the theoretically predicted rates as we increase the blending 
width K. 

We use two kinds of blending functions: a cubic spline 



B{x) 







-2x^ + 3x2 



1 



and a quintic spline 



B{x) 





6x^ 
1 



ISx-^ + lOx^ 



X < 0, 
< X < 1, 

X > 1, 

X < 0, 
0<x< 1, 
X > 1. 



(2.10) 



(2.11) 



We scale B{x) and B[x) and define the blending functions for the atomistic chains as 

Pr.= B(^ andk-=B(^] foTi = -N + l,...,N. 

Therefore, atoms with indices from —N + 1 to belong to the continuum region, from 1 to iC — 1 
belong to the blending region, and from K to N belong to the atomistic region. We note that B{x) 



has three bounded derivatives and hence it satisfies (2.5), whereas for B{x) the second derivative 



has a jump, hence the third derivative does not exist. Therefore, we expect that only (3 will yield 
the asymptotically optimal stability estimates for the B-QCF method (see [22j). 
For our interaction potential, we use the Morse potential 

l2 



(r) = [1 — exp(— a(r — 1))]^ 



(2.12) 



and we cut-off the interactions beyond the second nearest neighbor interactions. 

We apply a uniform expansion to the atomistic chain: yp := FeZ with Dirichlet boundary 
condition: 

u-AT+i = uj\[ = 0. (2-13) 

We then compute the critical strains of the atomistic and B-QCF models with different blending 
size K and fixed A^. The critical strains are defined as 



7™ := max{F > : L^iyc) is positive definite for all G G [1, F)} 
where w S {a, c, bqcf} denotes the respective model. 



(2.14) 



Remark 2.1. The stability bounds in Theorem 2.1 hold also for displacements u satisfying a ho- 
mogeneous Dirichlet boundary condition. To establish this, we note (1) that the bounds hold for 



constant displacements as well, and (2) that any function satisfying (|2.13) can be extended to a 



periodic function (possibly with a nonzero mean). Hence, Corollary |2 . 1| also holds for displacements 



u with homogeneous Dirichlet boundary conditions (2.13). 
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The computational results are shown in Figure [!} In Figure 1(a) we plot the dependence of the 
errors of quintic blending on K for different values of a. We see that the graph of th e er ror for the 
quintic blending is very close to the lower bound K~^''^ as given by (2.^ 



in Theorem 



2.1 Also, the 



error is lower for larger a, which is also in accordance with the theoretical results. Indeed, when a 
is large, the strength of the next-nearest neighbor interaction, 
neighbor interaction 



'2P, is small relative to the nearest 



bp, which contributes to a better stability of B-QCF according to (2.8). 



Figure 1 (b) shows the results of comparison of the cubic and the quintic blending. We see that 
the cubic blending produces the error that seems to decay slower, like K~^. On the other hand, 
the quantitative difference between cubic and quintic is not large on the example considered. To 
observe a significantly higher accuracy of the quintic spline, the computational domain size N has 
to be much larger. In addition, for larger a, N has to be even larger for the quintic blending to 
have advantage over the cubic blending. 



Absolute strain error for 1 D uniform expansion: N=2048 




iog,„(K) 
(a) Quintic spline blending 



Abs strain error for 1 D uniform expansion: N=4096,a=3 

—©—cubic spline 
—^quintic spline 




iog,„(K) 
(b) Quintic v.s. Cubic 



Figure 1. (a) The absolute critical strain errors for a ID uniform expansion. 
We set A^ = 1024, A7 = 1/A^^ where A7 is the strain increment used for testing 
stability, and 7^ and 7^*5"^^ are the critical strains for the atomistic and B-QCF 
models, respectively. The dashed line corresponds to the theoretical asymptote, 
(b) The absolute critical strain errors of quintic and cubic blending functions with 
N = 4096 and a = 3. The solid line corresponds to the theoretical asymptote. 



3. The B-QCF Operator in 2D. 

3.1. The Triangular Lattice. For some integer iV G N and e := 1/A^, we define the scaled 2D 
triangular lattice L to be 

1 1/2 
\/3/2 

where Cj, i = 1,2 are the scaled lattice vectors. Throughout our analysis, we use the following 
definition of the periodic reference cell 

n:=AQ{-N/2,N/2f and £ := L n O. 



L 



AgZ^, where Ag := [01,02] := e 
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We furthermore set 03 = (— l/2e, ^/?)/2eY , then the set of nearest-neighbor directions is given by 

Ni := {±ai,±a2,±a3}. 

The set of next nearest-neighbor directions is given by 

7V2 := {±6i,ib62,±&3}, where 61:= 01 + 02, b2:=a2 + as, and 63 = 03—01. 

We use the notation M := A/i U M2 to denote the directions of the neighboring bonds in the 
interaction range of each atom (see Figure [2]). 

We identify all lattice functions v : L — t- M^ with their continuous, piecewise affine interpolants 
with respect to the canonical triangulation T of M'^ with nodes L. 



000000 



o o o , • o o o 




o o 



o o o • o o o 



000000 
(a) Neighbor set 



(b) Domain decomposition 



Figure 2. (a) The 12 neighboring bonds of each atom, (b) The periodic reference 
cell £ := L n fi, the atomistic region f^a := Hex(ei?a), and the blending region 
fib := Hex(ei?fe) \ Qa- Here, N = 32, Ra = 3, Rb = 7, and K = A. 



3.2. The Atomistic, Continuum, and Blending Regions. Let Hex(i?) denote the closed 
hexagon centered at the origin, with sides aligned with the lattice directions 01,02,03, and di- 
ameter 2R. For Ra < Rb < iV € N, we define the atomistic, blending, and continuum regions, 
respectively, as 

fifl := Hex(ei?a), ^6 := Hex(ei?;,) \ fi^, and fie := clos (fi \ (fia U fib)) . 

We denote the blending width by K := Rb — Ra- Moreover, we define the corresponding lattice 
sites 

/:":=£nfia, £^=£nfife, and £'=:=£nfic. 

For simplicity, we will again use C as the finite element nodes, that is, every atom is a repatom. 
For a map v : L — t- M^ and bond directions r,s ^ J\f, we define the finite difference operators 

v{x + r)-v{x) r, r, f \ Dsv{x + r) - Dsv{x) 
DrV[x) := and DrL>sV[x) := . 

We define the space of all admissible displacements, U, as all discrete functions L — )• M^ which 
are fi-periodic and satisfy the mean zero condition on the computational domain: 

d2 



u 



-{- 



L 



: u{x) is fi-periodic and '^x&c''^(^) "*-*(■ 
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For a given matrix B G M^^^, det(B) > 0, we admit deformations y from the space 

3^5 := {y : L -> R^ : y{x) = Bx + u{x)\/x € L, for some u^U). 

For a displacement u € ZY and its discrete directional derivatives, we employ the weighted discrete 
l1 and l'^ norms given by 

/ xl/2 

||u||^2 := e ^^|n(j;)| , ||u||^c>o :=max|u(x)|, and 

\ xec J 

1/2 

\Dn\\(2 



\ x&C i=l J 



The inner product associated with l^ is 

(u, w) := e^ 2, u{x) ■ w{3 



xec 

3.3. The B-QCF operator. The total scaled atomistic energy for a periodic computational cell 
0, is 

^^(y) =f E E ^(Drvix)) = e^ E E [^(DaM^)) + HDbM^))] , (3.1) 

xeCreAf xeC i=l 

where (j) £ C'^{M.'^), for the sake of simplicity. Typically, one assumes (/)(r) = <^(|r|); the more general 
form we use gives rise to a simplified notation; see also 35 . We define (p'{r) £ M^ and ^"(r) E M^^^ 
to be, respectively, the gradient and hessian of (p. 

The equilibrium equations are given by the force balance at each atom, 

F^{x;y) + f{x;y) = 0, for x G C, (3.2) 

where /(x; y) are the external forces and F"'{x; y) are the atomistic forces (per unit area e^) 

1 d8-{y) 



F%x;y):-- 



e^ dy{x) 

3 , 3 



i=l i=l 

Again, since u = y — y^, where ysix) = Bx, is assumed to be small, we linearize the atomistic 



equilibrium equation (3.2) about y^: 

(L^u^)(x) = /(a;), for x G £, 
where (L^u) (x), for a displacement u, is given by 

3 3 

(L^u) (x) = -Y,4>"{Ba^)Da,DaM^ " a^) -J2^"{Bbi)Db,DbM^ - h), for x e £. 
We use the Cauchy-Born extrapolation rule to approximate the nonlocal atomistic model by a 



local continuum Cauchy-Born model 30,41,43]. Using the bond density lemma ^35^ Lemma 3.2] 
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(see also [40]), we can write the total QCL energy (the discretized Cauchy-Born energy) as a sum 
of the bond density integrals 

where dry{x) := ^y{x + tr)\t=o denotes the directional derivative. We compute the continuum 

force 

1 dS'^ 

and linearize the force equation about the uniform deformation y^ to obtain 

(L=u^)(x) = /(a;), for x G £. 

To formulate the B-QCF method, we let the blending function /3(s) : E? — )• [0, 1] be a "smooth", 
ri-periodic function. Then, the (nonlinear) B-QCF forces are given through a convex combination 

of F^[x; y) and F^(x; y): 

F'^^'^ix^y) := P{x)F^{x;y) + (1 - /3(x))F^(x;y), 
and linearizing the equilibrium equation F^ + / = about ys yields 

(^bqcf^bqcf)^^) = /(x), for X G £, 

(3.4) 
where (L^^^M(x) = /3(x)(L^u)(x) + (1 - /3(x))(L=u)(x). 

The 2D blending function in our computational experiments will be defined radially using cubic 
and quintic splines: 

f^i-) ■■= B ( ^S^) and ^(x) := 5 ^ ^^^ - """ 



eRb - eRa J \ eRb - eRa 



where B{x) is given by (2.10) and B{x) is given by (2.11). The function j3{x) has the smoothness 



and satisfies 2D versions of the scaling bounds (2.5) needed for Theorem 3.1 below, whereas (3{x 



does not have a bounded third derivative. We therefore can expect that /3(x) will give a larger 
error asymptotically as compared to (3{x). 

3.4. Positivity of the B-QCF operator in 2D. Necessary and sufficient conditions for L^^"^^ to 
be positive-definite are given in 122J. To make this paper more concise, we only state the conclusions 
without proof. First, we state an upper bound for {L^^'^^u,u): 

Theorem 3.1. Suppose that /3 G C^ and satisfies the scaling bounds (2.5); then, 

where 

7bqcf := 7 - C [K-'/^rI^'\ log{Rb/N)\'/^] , (3.5) 

where C is a generic constant independent of N, and 7 is the coercivity constant for the operator 
L: 

3 

(Lu,u) := (L"u,u)-e^^ J]/3(x-a2)|Z)a.^a,+iU(x-ai-a2)|^. > 7||Z)u||| Vu G Z^. 

One can see very clearly that, whenever A^ is polynomial in Rb and K ^ RJ , then L^'^ can be 
expected to be coercive. Both are natural and easy to achieve. We can thus deduce the following 
result for the coercivity of L '^" : 
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Corollary 3.1. Suppose that L is positive- definite and that the blending function /3 G C^ and 
satisfies the scaling bounds (2.5). Let the number of atoms Ra along the radius be of order N'^ with 



< a < 1 . If the blending width K satisfies 

' logiV|l/^ 



iv:> < 



a 



|logAr|i/57v°/5, 
A^i/5 a = 1, 



0, 
< a < 1, 



then the B-QCF operator L '^'^ is positive- definite. 
Remark 3.1. 



(a) The stability result of Theorem 3.1 , and hence of Corollary 3.1 , is conditional on the stability 
of the operator L. In 22 we show that L is indeed stable whenever nearest-neighbor 
interactions dominate. 

Moreover, based on the analysis and numerical experiments in [35i for a similar linearized 
operator, we expect that the region of stability for L is the same as for L" as N^ Ra, Rb — )• oo. 



We therefore expect that the result of Theorem 3.1 holds (up to a controllable error) if 



coercivity of L is replaced by coercivity of L^ in the hypothesis. 



(b) One can apply the argument of Remark 2.1 to conclude that the results of Theorem 3.1 



and Corollary 3.1 are valid for homogeneous Dirichlet boundary conditions as well. 



By constructing a radial counterexample similar to our ID counterexample, we can observe that 
our conditions in Corollary |3 . 1 1 are essentially necessary. 

Theorem 3.2. Suppose that L^ is positive- definite and that the blending function /3 G C^ and 
satisfies the scaling bounds (2.5). The number of atoms Ra along the radius is of order N" with 
< a < 1. If the blending width K is K <^ N'^'^, then the B-QCF operator L^^'^^ cannot be 
positive-definite and we can construct a radial counterexample in this case. 

We note that there is a gap between the necessary and sufficient conditions for < a < 1. In 
addition, we have no necessary condition for a = 0, which corresponds to a fixed atomistic core 
independent of the reference cell fi. 

3.5. 2D numerical experiments for B-QCF operators. In this subsection, we will continue 
the numerical experiments for the 2D B-QCF models to verify the theoretical findings by comparing 
the decay rates of the error in critical strain as computed by B-QCF with the theoretically predicted 
rates as we increase the blending width K. 

(1) Uniform expansion. 

We first consider the simplest 2D deformation: we apply a uniform expansion y{x) = Ex 
with 

to the perfect lattice C with Dirichlet boundary condition: 

u{x) = Vx e d^. (3.6) 

Then we compute the critical strains 7 of the atomistic and B-QCF models with different 
blending region width K. 

We note that the 2D conclusions also depend on the size of the atomistic region. Therefore 
we let Ra = K^'^ in order to narrow the dependence only to the blending width K. Then 
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the asymptotical term in (|3.5|) for sufficiently large A'^ is approximately 



which means that the error in 7bqcf is systematically improvable. 

The choice of scaling Ra = K^'^ is motivated by the results in (iTJ which indicate that, 
generically, one should expect an 0{R~^) error in the regions of stability between the infinite 
lattice atomistic model and the atomistic model in a domain with radius Ra- 

The critical strains are defined as 

7™ := max {7 > : V"{ByL) is positive definite for 7 G [0,7)}, (3-7) 

where w G {a, bqcf} denote the models. Here we use the MATLAB function eigs [29] to 



compute the smallest eigenvalue of the symmetric part of L" (i?x) and thus determine the 
positive-definiteness of L^(Bx). 

We also define the increment of the strain 7 in each step by A7. The results in |10| 
[l4||T7||35] suggest that the theoretical increments be of order 0{N~'^) (at least, for finding 
the critical strain of a uniform lattice), and we set A7 = 10^^ which is sufficiently small 
considering N = 200 or 300 in our experiments. 



Abs strain error for 2D uniform expansion: N=300, R =K ' 




Abs strain error for 2D uniform expansion: N=300, R =K ,a=3 




(a) Quintic spline blending 



(b) Quintic v.s. Cubic 



Figure 3. (a) The absolute critical strain errors for the 2D uniform expansion. 
We set A'^ = 300, and we denote the critical strains for the atomistic and B-QCF 
models by 7^^ and 7"'^^ , respectively. The dashed line corresponds to the theoretical 
asymptote, (b) The absolute critical strain errors for the quintic and cubic blending 
functions with A^ = 300 and a = 3. The solid line corresponds to the theoretical 
asymptote. 

We plot the difference of the critical strains with different blending width K in Figure |3j 
The numerical critical strain errors in the left figure approach the analytical asymptote as 

K increases. There are larger fluctuations of errors as compared to the ID case, which is 

3/5 
likely due to round-off errors in calculating K = Ra ■ Thus, the slopes of the errors with 



quintic blending agree with the theoretical prediction in Theorem 3.1 Also, similarly to 
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the ID results, the error is smaller when the nearest neighbor interaction dominate (that 
is, when a is large). 

Although the slope of the errors with cubic blending seems to be one half order less 



than that with quintic blending (see Figure |3(b) ), the computed errors for cubic blending 
are slightly smaller for the relatively small value of A^ considered. We expect that for a 
sufficiently large system, the quintic blending would be more accurate. In addition, the 2D 
errors for uniform expansion are similar to the ID results. This is reasonable since the 2D 
uniform expansion is similar to the ID deformation. 
(2) Uniform shear deformation. 

We now investigate stability of B-QCF under shear deformation. We apply a y-directional 



shear deformation to the hexagonal lattice Vt with Dirichlet boundary conditions (3.6). The 
y-directional shear is y{x) = Bx with 



B 



The critical strain errors between the B-QCF and atomistic models with the quintic blending 
are plotted in Figure |4| 
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(b) Error vs Ra 



Figure 4. The relative critical strain error for the y-directional shear deformation. 
7^^ and 7'^'^'^^ are the critical strains for the atomistic and B-QCF models respectively. 
The dashed line corresponds to the theoretical asymptote. The fluctuations in the 
plotted error for N = const seems to be due to round-off errors in calculating Ra 
and K. 

In FigureH^we plot the critical strain errors in the following three regimes: (1) A increases, 
Ra = const, K = const, (2) all three parameters increase, and (3) A = const, Ra and K 
increases. The results indicate that the error in this case depends on A, but does not depend 
on Ra or K. This means that, for shear deformations, the local continuum approximation 
and its finite element coarse-graining contributes most of the error. 

We explain such a qualitative difference between the uniform expansion and the shear 
deformation in the following way. For the uniform expansion the onset of instability is 
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(3) 



due to competition of interaction of the nearest neighbors (NNs), contributing to stabihty, 
and the second nearest neighbors (NNNs), contributing to instabiUty. On the other hand, 
for the shear deformation the onset of instability is primarily due to competition between 
elongated and compressed NN bonds. Therefore, for the uniform expansion it is important 
to reduce the interface error which distorts the NNN interaction, whereas in shear defor- 
mation the NNN interactions do not contribute significantly to stability errors. Since, for 
NN interaction, the atomistic, Cauchy-Born and B-QCF models are identical, the stability 
error only depends on the domain size. 
Regions of stability. 

We now combine the uniform expansion and shear deformation together and study the 
stability region of L^^^^ for a general class of homogeneous deformations. We consider the 
following group of deformations which involve general shear, expansion, and compression. 



B 



1 + s 




0.1 

1 + r 



Applying these specific homogeneous deformations to the hexagonal lattice in the reference 
cell and again using the Dirichlet boundary condition, we plot the stability regions (regions 
where the operators are positive definite) in Figures] 
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Figure 5. The stability regions of the different models. These closed curves are the 
boundaries of the stability regions for the atomistic, B-QCF, and the local continuum 
models, respectively. The curves with indicators are for QCF. 



We observe that the stability regions of the B-QCF model with different blending sizes 
are all proper subsets of the atomistic model. In addition, the fully atomistic and continuum 
models are very close to each other, which agrees with the stability analysis of the perfect 



lattice 17 . Also, when a increases, which means the next-nearest neighbor interactions 



become less important, the difference becomes smaller. 
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There is a visible difference in the stability regions between the QCF model and the exact 
atomistic model, whereas the difference between the B-QCF model and the atomistic model 
is almost not seen. This implies that using a blending region can significantly improve the 
stability properties of the approximation models. 
(4) Stability of micro-cracks. 

The experiments that we have reported up to this point were based on perfect lattices. 
Now we apply the B-QCF model to lattices with local defects. 

The atomistic system is as follows. There is a micro-crack in the center of the domain 
Q, with length 5, i.e., 5 atoms are removed from the lattice (see Figure pi). We impose a 

vertical stretching B = ^ ^ , on the lattice and compute the critical strains i'^ > 

1 + 7 

beyond which the system loses stability. 



Figure 6. The stable equilibrium configuration of the micro-crack with crack 
length= 5 and 7 = 0.001, and the i'^ norm of the force residual is of order O(10~^^). 

We computed the critical strain 7^^ in the following way. Given 7 > 0, we use Newton's 
iteration method to solve the following force equations for y"^ with the initial guess yp = Bx: 

F"'(x;yT)=0 for xen\dn. 

We set the tolerance for the i'^ norm of the force residual of the Newton's iteration to 
be 10~^. To prevent the configuration from "jumping out" of the local energy well corre- 
sponding to the defect under consideration, we require at each step the i"^ norm of force 
residual to be less than 100 and the positive-definiteness of -L"'(y), where y is the current 
configuration. If any of the two requirements is not met, then the current 7 is regarded as 
an unstable strain. When the force residual is smaller than the tolerance, the configuration 
y* is thought to be in its equilibrium of the local energy well. Then we check the positive 
definiteness of corresponding operator L"(y*) with the equilibrium configuration y*. The 
nonlinear critical strain is thus defined as 

7^^ := max{7 > : L""(y*) is positive definite for 7 G [0, 7)}. 

The plot of critical strain for the B-QCF models are shown in Figure [7j We observe 
the nonlinear error decays much faster than the theoretical predicted rates and it can 



reach the strain increment A7 = 10~ . This phenomenon has been observed in 35 and 
is likely related to superconvergence of local quantities of interest. The indicator of the 
superconvergence is the concentration of the critical eigenmode corresponding to 7*^ near 
the defect, which is illustrated in Figure |8| 
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Figure 7. The nonlinear critical strain error for vertical stretching. We set A^ = 
200, crack length= 5, and Ra = TaaK.{K^,6}. 7*^, 7^'^^^ are the critical strains for 
the atomistic and B-QCF models, respectively. The dashed line corresponds to the 
theoretical asymptote. 



■'■/I' \', 



Figure 8. The zoomed-in critical eigenvector of critical strain of vertically stretch- 
ing a micro-crack. We set N = 200, crack length= 5, strain increment A7 = 10"^^ 
and a = 4.2. 



We also study the relative errors of the critical strains for two different choices of the 
blending width, K = 2 and K ~ RJ + 2. Motivated by the analysis in '35), the size of 
the atomistic core is chosen to be Ra = vN . According to Figm'e k)l the relative errors for 
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3/5 

K ~ Ra + 2 are approximately 10 times smaller than those for K = 2. But both graphs 
decay rapidly as N increases. The rate of decay appears to be quadratic. 
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Figure 9. The errors of the critical strains of vertically stretching a micro-crack 
in logj^o scale plot. We set crack length= 5, and a = 4.2. 7^, 7^'^^^ are the critical 
strains for the atomistic and B-QCF models, respectively. The dashed line is the 
theoretical asymptote. 



4. The Accuracy of B-QCF 

In the previous sections, we investigated the positivity of the B-QCF operator. One motivation 
for this study is that these experiments fill a gap in our error analysis of the B-QCF method |20|. We 
now briefly review these results and then include some numerical experiments demonstrating the 
superior accuracy of B-QCF over other a/c coupling schemes that we have investigated previously 
in [28]. 

4.1. Implementation of the B-QCF method. Let V C L be a set of vacancy sites and Ly := 
L \ V the corresponding lattice with defects. Let B € M?^'^ be the applied far- field strain. We 
consider the atomistic problem 

y"" e argmin{£:''(y) : y : Ly ^ ^"^jVix) ~ Bx as |^| -^ 00}. (4.1) 

We remark that one must carefully renormalize £^ in order to rigorously make sense of this problem; 
for the details. 



20 



see e.g. 

We wish to approximate this problem with a practical variant of the B-QCF method. To that 
end, we choose Ra,Rb = Ra + K,N E N in such a way that all vacancy sites are contained in 
the atomistic region Oq, which is a hexagon with side length Ra- The blending region is defined 
analogously. The full computational domain is given by fi, which is a hexagon with side length A^. 
We triangulate il in such a way that it matches the canonical triangulation of the triangular lattice 
in Q. 
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Let Th denote the set of triangles, let Af^ denote the nodes of the triangulation, and let AfJ^^*^ '■= 
Mh \ (V U (90) denote the free nodes. 

Let P^ denote the space of all functions f/j : J7 — )• M^, that are continuous and piecewise afiine 
with respect to the triangulation Th- The space of admissible trial functions is then given by 

Yh := [Vh G Pft : Vhix) = Bx for x^dVL). 

Each deformation yh G Y/^ is understood to be extended by Bx outside of Q. and thereby gives rise 
to an admissible atomistic configuration. 

We define the discretized Cauchy-Born energy functional as 

£'{yh)-= ^ vol(T)VF(Vy^|T), 

TdTh 

where vol(r) in 2D is the area of the triangle T. We can define the discretized B-QCF operator, 
for a given blending function /3, as follows: 

F'^^^^x; Vh) := (1 - /3(x))^^ + /3(x)f^ for x G <-. 



dy{x) 



y=yh 



dzh{x) 



Zh=yh 



In the B-QCF method, we aim to find a solution y^^^^ G Y^ satisfying 

F^'i'^ix; y^'^'^) = Vx G AA|■'^^ (4.2) 

We remark that this method has essentially five approximation parameters that must be chosen 
carefully: the atomistic region size Ra, the blending width K, the computational domain size N, 
the blending function f3, and the finite element mesh Th- 



4.1.1. Practical considerations. To implement (4.2) in practice, we need to specify further details 
of the method: 

(1) In our choice of blending function, we deviate from the optimal choice of a C^'^-blending 

function and instead choose only a C^'^ blending function, which is more easily constructed. 

To be precise, we choose the blending function proposed in 28 , which minimizes ||V^/3||j;^2, 



or a discrete variant thereof, in a precomputation step (see [28] for the details). We have 
seen in Figure [T] and Figure |3] that using the less regular blending function makes only 
a small difference for the stability results when the size of the computational domain is 
moderately large, and it guarantees that all our results discussed in section |4.2| are still 
valid. 
(2) In addition to the blending region Clf, we ensure that two additional "layers" of atoms outside 
of it belong to J\fh- This makes the implementation of the atomistic force contribution in 



(4.2) straightforward. 

Moreover, we ensure that the vacancy sites do not affect the forces on atoms x where 



f3{x) y^ 0. This ensures that all the Cauchy-Born force contributions in (4.2) are the correct 
Cauchy-Born forces. 
(3) To obtain an appropriate initial guess for the B-QCF solutions, we first solve the corre- 
sponding energy-based blended QCE method (B-QCE) [28] with the same approximation 
parameters, using a preconditioned line search method. The details are described in [28|. 
The B-QCE solution is then taken as a starting guess for the B-QCF Newton iteration to 



solve ( |4.2D . 

We remark, that the Jacobian matrix of the B-QCF operator is straightforward to as- 
semble from the Hessians of the atomistic and Cauchy-Born energy. Nevertheless, for large 
3D simulations, more sophisticated solution methods may be required. 
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(4) We are now only left to choose the remaining approximation parameters Ra,K,N and the 
mesh Th- 



4.2. Error versus computational cost. We briefly review the main ideas of our analysis in |20| 
without technical details. A first key result is that if the atomistic solution is stable {5'^£^{y'^) is 
positive definite) and the linearized B-QCF operator 5F^'^'^^{-; Bx) is positive definite, then choosing 
Ra and K sufficiently large implies that 5F^'^"{-]y^) is also positive definite, that is, the B-QCF 
method is stable under these conditions. To achieve this in practice, we need to choose K^ ^ Ra 



(recall from section 4.1.1 that we have chosen a sub-optimal /3). 

From this stability result, we can deduce the existence of a B-QCF solution in a neighborhood 
of the atomistic solution, and an error estimate in terms of the best approximation error (the 
best approximation of y^ from the finite element space Y^). and of the modeling error (the force 
discrepancy of the B-QCF and atomistic models). We estimate the error in the strain Vy^ — Vy^^'^ 
in terms of the "smoothness" of y^, which is measured in terms of bounds on the derivatives V^y^. 
The derivatives of the discrete functions y^ are understood as derivatives of a smooth interpolant. 



(See 20 for the details. 



Dropping an unimportant term for the sake of readability, our error estimate reads 

llVy^ - Vy^'^^IU^d,^) < C^*'^^(c'^||v3y-||^.(ij2\,^) + || W2y-|U.(^\,^) + ||Vy'^||i2(M2\^)) , (4.3) 

where ||V^y'^||i2(R2\^^) measures the modeling error, ||/iV^y^||j^2(f7\(^^) the finite element discretiza- 
tion error and ||Vy^||/^2(ig2\^-) the error in the far-field due to the artificial boundary condition (the 
two latter errors comprise the best approximation error). The domains Wa,w are slightly smaller 
hexagonal subsets of, respectively, f^a and f], with comparable side lengths. 

In addition, C^^'"^^ is a stability constant that is uniformly bounded for Ra <^ K^, and 

C^ ■.= K-^/^Rl/^log\Ra/N\ 

is a /3-dependent prefactor, which arises from a crucial inequality, ||V(/3t;)||^2 < C^||Vt;||j;^2, in the 
consistency analysis of B-QCF. 

We choose K ^ Ra and N a polynomial of Ra (we will see momentarily why this is natural), 
then C^ is uniformly bounded and in addition, we choose Ra ^ K^ , which we require for stability. 
With this choice, it is easy to see that C^ \\\/^y''^\\ i2(^-^2\i^^-^ < \\hV'^y^\\i^2^^Q^\^^-^ (recall that we are 
working in units where atomic spacing is 1), and hence we can simply ignore the modeling error 
term from now on. 

We recall from [28| that the atomistic method (ATM) is given by the B-QCF method with /3 = 0. 
We also recall the corresponding error estimates (dropping less important terms) for the atomistic 
(ATM) and the B-QCE methods [20l[28] 



||Vy- - Vy-*-|U2(M2) < \\Vy-\\LHR^\^), (4.4) 

||Vy" - Vy^^"li2(M2) < ||v2/3||i2(M2\^j + \\hV^yl mn\u.^) + \Ny"\\LHR^\u.)- (4.5) 

To better understand the best approximation error, we need to understand the regularity of y^. 
Since the problems only involve defects with zero Burgers vector, it is reasonable to assume based 
on linear elasticity, that 

|V^y'^(x)| ~ \x\-^-\ 

Having this explicit knowledge about the elastic field, we can optimize our choice of finite element 
triangulation. Using the construction in J35j and also used successfully in our B-QCE experiments 
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in 28 , we obtain a triangulation Tii (as a function of i?5 and N)^ for which the following estimate 
holds: 

Wh^h^WL^iSA^^) + l|Vy"||L2(K2\^) < R-^ + N-\ 

Thus, we choose A^ ~ i?^ to balance these two error contributions. 

Finally, we note that, with this construction, the number of degrees of freedom in Y/^, DoF := 
dimY/i = 2^J\fl^^^ is approximately equal to DoF ~ R^. (In particular, the number of degrees of 
freedom in the atomistic, blending and continuum regions are comparable.) 

In summary, choosing K ~ Ra, N ~ i?^, the blending function /3 according to the construction 
proposed in |28|, and the finite element mesh according to the construction proposed in 35 , we 



obtain from (4.3) the error estimate 

l|Vy--Vy^'=fb2(K2)<DoF-i. (4.6) 



We note that N = Ra in the ATM method, and consequently we obtain from (4.5) 

||Vy--Vy-*-|U2(M2)<DoF-V2; (4.7) 

thus demonstrating an improved rate of convergence for the B-QCF method in comparison with 
the ATM method. 

We remark that this is optimal for Pl-finite element type coarse-graining schemes, as the mod- 
eling error is in fact dominated by the finite element error. In particular, it is a substantial 
improvement over the B-QCE method, for which the corresponding error estimate obtained from 



( 4.4D is 

We note that the B-QCE method can be shown to have a higher rate of convergence than the ATM 
method for defects with nonzero Burgers vector (such as dislocations) which have a lower rate of 
decay. The finite element coarse-graining of the B-QCE method can more efficiently approximate 



the larger region where the strain gradient is significant; 20 281 for the details. 



4.3. Numerical rates. We test our analytical predictions against the two numerical examples, 
for which we already tested the B-QCE method in [28] . In both examples, we choose the Morse 
interaction potential 

</,(r) = e-2"('-i)-2e-"('-^), 
with stiffness parameter a = 4. 

We compare the B-QCF method with a pure atomistic computation on a finite domain, with the 
QCE and B-QCE methods (cf. [28] for a detailed description of these three methods) and with the 
pure QCF method, which is simply the B-QCF method with K < 1 (i.e., f3{x) £ {0, 1}). 

Finally, we have also included a highly optimized B-QCE variant where we choose K ~ R^ 
and N ~ i?^, which is a very unexpected scaling, but yields improved errors in the preasymptotic 



regime; see 28, Remark 4.3]. We denote this method by B-QCE+ in the error graphs. 



4.3.1. The di-vacancy example. We choose the vacancy set V = {0, ei} and the macroscopic strain 

/1.03 0.3 \ 
^"VO.O 1.03J"^°' 

where Bq is a minimizer of W (3% uniform stretch and 3% shear from ground state). The setup of 



the B-QCF method for the di-vacancy problem is shown in Figure 10 



In Figure 11, we plot the degrees of freedom (DoF) against the error in the energy- norm for 
the various a/c coupling methods that we consider. As predicted by our analysis, the B-QCF 
method clearly outperforms all other methods, with the exception of the QCF method, which is 
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Figure 10. Setup of the B-QCF method for the di-vacancy example, for a specific 
choice of approximation parameters, shown in deformed equihbrium. The size/color 
of the atoms in the center correspond to decreasing values of (1 — /3{x)). 

Convergence rates for the di-vacancy example 
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Figure 11. Plots of computational cost (DoF) versus error in the energy- norm 
for various a/c coupling methods approximating the di-vacancy problem described 
in section 14.3.11 



barely distinguishable from the B-QCF method in this graph. Unfortunately, we cannot offer a 
satisfactory theory for the QCF method at present. 

We also remark that, due to the high consistency error committed in the interface region, the 
B-QCE does not even outperform a plain atomistic computation in this particular example. (But 
it will clearly outperform the fully atomistic method (ATM) in the micro-crack example, where the 
elastic field is much more significant.) 
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Figure 12. Setup of the B-QCF method for the micro-crack example, for a specific 
choice of approximation parameters, shown in deformed equihbrium. The size/color 
of the atoms in the center correspond to decreasing values of (1 — /3{x)). 
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Convergence rates for the microcrack example 
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Figure 13. Plots of computational cost (DoF) versus error in the energy- norm 
for various a/c coupling methods approximating the micro-crack problem described 
in section 14.3.21 



4.3.2. The micro-crack example. In the micro-crack example, we choose the vacancy set V = 
{— 5ei, . . . , 5ei} and the macroscopic strain 

_ A.O 0.03\ 
^ - VO.O I.03J ■ ^0' 

where Bq is a minimizer of W (3% tensile stretch and 3% shear from ground state). The setup of 



the B-QCF method for the micro-crack problem is shown in Figure 12 
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In Figure 13 we plot the degrees of freedom (DoF) against the error in the energy- norm, for the 
various a/c couphng methods that we consider. In this example the picture is less clear than in the 
di-vacancy example due to a more significant preasymptotic regime, which is caused by the more 
significant deformation admitted by the microcrack. In the preasymptotic regime we observe that 
the QCE and B-QCE methods perform much better than expected, but eventually fall back to the 
predicted rates. By contrast, the B-QCF and QCF methods display clear systematic convergence 
at the predicted rate throughout. 

We also note that, in this example, the B-QCE+ method performs comparable to the B-QCF 
and QCF methods, at least in the preasymptotic regime accessible in the experiment. 

5. Conclusion 

We have formulated an atomistic-to-continuum force-based coupling, which we call the blended 
force-based quasicontinuum (B-QCF) method. In this paper, we numerically studied the stability 
as well as accuracy of the B-QCF method. We computed the critical strain errors between the 
atomistic and B-QCF models with different sizes of the blending region under different types of 
deformations. 

The main theoretical conclusion in I22i is that the required blending width to ensure coercivity 
of the linearized B-QCF operator is surprisingly small. For both ID and 2D uniform expansion, 
the computational results of the linearized operators perfectly match the analytic predictions. In 
addition, the stability for a general class of homogeneous deformations of the 2D B-QCF operator 
becomes almost the same as that of the atomistic model by using a very small blending region, 
in contrast to the fact that the stability region of the force-based quasicontinuum (QCF) method, 
that is, the B-QCF method without blending region, is just a proper subset of the fully atomistic 
model. However, the critical strain error for the B-QCF operator applied to shear deformation 
seems to only linearly depend on the system size and is thus insensitive to blending width. 

For the problem of a microcrack in a two-dimensional crystal, we studied the nonlinear stability 
of the B-QCF operators. The critical strain error decays faster than the prediction, and it can be 
as small as the strain increment. However, we find that the error increases a little bit when the 
blending size becomes larger, which is possibly due to round-off error. 

Moreover, we implemented a practical version of the B-QCF method. We briefly reviewed the 
accuracy results in terms of computational cost [20]. The numerical experiments, di-vacancy and 
microcrack demonstrate the superior accuracy of B-QCF over other a/c coupling schemes that we 
have investigated previously in [28j . 

The BQCF method with a surprisingly small blending region is an appealing choice for numerical 
simulations of atomistic multi-scale problems as it is always consistent and can be guaranteed by 
both theory and benchmark testing to be positive definite when the fully atomistic operator is 
positive definite. 
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